{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using VMLS\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chapter 6\n",
    "# Matrices\n",
    "### 6.1 Matrices\n",
    "**Creating matrices from the entries.** Matrices are represented in Julia as $2$-dimensional arrays. These are constructed by giving the elements in each row, separated by space, with the rows separated by semicolons. For example, the $3 × 4$ matrix\n",
    "$$\n",
    "A = \\begin{bmatrix}\n",
    " 0 & 1 & −2.3 & 0.1\\\\\n",
    " 1.3 & 4.0 & −0.1 & 0\\\\\n",
    " 4.1 & −1.0 & 0 & 1.7\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "is constructed in Julia as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Float64,2}:\n",
       " 0.0   1.0  -2.3  0.1\n",
       " 1.3   4.0  -0.1  0.0\n",
       " 4.1  -1.0   0.0  1.7"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [0.0 1.0 -2.3 0.1;\n",
    "    1.3 4.0 -0.1 0.0;\n",
    "    4.1 -1.0 0.0 1.7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Here, `Array{Float64,2}` above the array tells us that the array is $2$-dimensional, and its entries are `64-bit` floating-point numbers.) In this example, we put the different rows of the matrix on different lines, which makes the code more readable, but there is no need to do this; we get the same matrix with"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Float64,2}:\n",
       " 0.0   1.0  -2.3  0.1\n",
       " 1.3   4.0  -0.1  0.0\n",
       " 4.1  -1.0   0.0  1.7"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [0 1 -2.3 0.1; 1.3 4 -0.1 0; 4.1 -1 0 1.7]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Julia function `size(A)` gives the size, as a tuple. It can also be called as `size(A,1)` or `size(A,2)` to get only the number of rows or columns. As an example, we create a function that determines if a matrix is tall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m, n = size(A)\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(A,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(A,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tall(X) = size(X,1)>size(X,2);\n",
    "tall(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the function definition, the number of rows and the number of columns are combined using the relational operator `<`, which gives a $Boolean$.\n",
    "\n",
    "**Indexing entries.** We get the $i, j$ entry of a matrix `A` using `A[i,j]`. We can also assign a new value to an entry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.1"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[2,3] # Get 2,3 entry of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Float64,2}:\n",
       " 0.0   1.0   7.5  0.1\n",
       " 1.3   4.0  -0.1  0.0\n",
       " 4.1  -1.0   0.0  1.7"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1,3] = 7.5; # Set 1,3 entry of A to 7.5\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Single index indexing.** Julia allows you to access an entry of a matrix using only one index. To use this, you need to know that matrices in Julia are stored in column-major order. This means that a matrix can be considered as a one - dimensional array, with the first column stacked on top of the second, stacked on top of the third, and so on. For example, the elements of the matrix\n",
    "$$\n",
    "Z =\n",
    "\\begin{bmatrix}\n",
    "−1 & 0 & 2\\\\\n",
    "−1 & 2 & 3\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "are stored in the order\n",
    "$$−1, −1, 0, 2, 2, 3.$$\n",
    "With single index indexing, `Z[5]` is the fifth element in this sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = [ -1 0 2; -1 2 -3];\n",
    "Z[5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is very much not standard mathematical notation, and we would never use this in VMLS. But it can be handy in some cases when you are using Julia. \n",
    "\n",
    "**Equality of matrices.** `A == B` determines whether the matrices `A` and `B` are equal. The expression `A .== B` creates a matrix whose entries are $Boolean$, depending on whether the corresponding entries of `A` and `B` are the same. The expression `sum(A .== B)` gives the number of entries of `A` and `B` that are equal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "false"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = copy(A);\n",
    "B[2,2] = 0;\n",
    "\n",
    "A == B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 BitArray{2}:\n",
       " true   true  true  true\n",
       " true  false  true  true\n",
       " true   true  true  true"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A .== B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Row and column vectors.** In Julia, as in VMLS, $n$-vectors are the same as $n × 1$ matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1×3 Array{Float64,2}:\n",
       " -2.1  -3.0  0.0"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [ -2.1 -3 0 ] # A 3-row vector or 1x3 matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -2.1\n",
       " -3.0\n",
       "  0.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = [ -2.1; -3; 0 ] # A 3-vector or 3x1 matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output reveals a small subtlety that generally won’t affect you. You can see that b has type `Array{Float64,1}`. The final `1` means that it is a $1D$ array; `size(b)` gives `(3,)`, whereas you might think it would or should be `(3,1)`. So you might say that in Julia, $n × 1$ matrices are n-vectors. This is why we say above that n-vectors and $n × 1$ matrices are almost the same in Julia. \n",
    "\n",
    "**Slicing and submatrices.** Using colon notation you can extract a submatrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Int64,2}:\n",
       " -1   0   1  0\n",
       "  2  -3   0  1\n",
       "  0   4  -2  1"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ -1 0 1 0 ; 2 -3 0 1 ; 0 4 -2 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Int64,2}:\n",
       " 1  0\n",
       " 0  1"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[1:2,3:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is very similar to the mathematical notation in VMLS, where this submatrix would be denoted $A_{1:2,3:4}$. You can also assign a submatrix using slicing (index range) notation. A very useful shortcut is the index range `:` which refers to the whole index range for that index. This can be used to extract the rows and columns of a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Int64,1}:\n",
       "  1\n",
       "  0\n",
       " -2"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[:,3] # Third column of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       "  2\n",
       " -3\n",
       "  0\n",
       "  1"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A[2,:] # Second row of A, returned as column vector!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In mathematical (VMLS) notation, we say that `A[2,:]` returns the transpose of the second row of `A`. As with vectors, Julia’s slicing and selection is not limited to contiguous ranges of indexes. For example, we can reverse the order of the rows of a matrix `X` using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Int64,2}:\n",
       "  0   4  -2  1\n",
       "  2  -3   0  1\n",
       " -1   0   1  0"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = copy(A) #For example\n",
    "m = size(X,1)\n",
    "X[m:-1:1,:] # Matrix X with row order reversed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia’s single indexing for matrices can be used with index ranges or sets. For example if $X$ is an $m × n$ matrix, `X[:]` is a vector of size $mn$ that consists of the columns of `X` stacked on top of each other. The Julia function `reshape(X,(k,l))` gives a new $k × l$ matrix, with the entries taken in the column-major order from `X`. (We must have $mn = kl$, i.e., the original and reshaped matrix must have the same number of entries.) Neither of these is standard mathematical notation, but they can be useful in Julia."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Int64,2}:\n",
       " 1  -3\n",
       " 2   0\n",
       " 1  -2"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [ 1 -3 ; 2 0 ; 1 -2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Int64,1}:\n",
       "  1\n",
       "  2\n",
       "  1\n",
       " -3\n",
       "  0\n",
       " -2"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Int64,2}:\n",
       " 1   1   0\n",
       " 2  -3  -2"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reshape(B,(2,3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "ename": "DimensionMismatch",
     "evalue": "DimensionMismatch(\"new dimensions (3, 3) must be consistent with array size 6\")",
     "output_type": "error",
     "traceback": [
      "DimensionMismatch(\"new dimensions (3, 3) must be consistent with array size 6\")",
      "",
      "Stacktrace:",
      " [1] (::getfield(Base, Symbol(\"#throw_dmrsa#193\")))(::Tuple{Int64,Int64}, ::Int64) at ./reshapedarray.jl:41",
      " [2] reshape(::Array{Int64,2}, ::Tuple{Int64,Int64}) at ./reshapedarray.jl:45",
      " [3] top-level scope at In[35]:1"
     ]
    }
   ],
   "source": [
    "reshape(B,(3,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Block matrices.** Block matrices are constructed in Julia very much as in the standard mathematical notation in VMLS. You use `;` to stack matrices, and a `space` to do (horizontal) concatenation. We apply this to the example on page [109](https://web.stanford.edu/~boyd/vmls/vmls.pdf#equation.6.1.1) of VMLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Int64,2}:\n",
       " 0  2  3  -1\n",
       " 2  2  1   4\n",
       " 1  3  5   4"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [ 0 2 3 ]; # 1x3 matrix\n",
    "\n",
    "C = [ -1 ]; # 1x1 matrix\n",
    "\n",
    "D = [ 2 2 1 ; 1 3 5]; # 2x3 matrix\n",
    "\n",
    "E = [4 ; 4 ]; # 2x1 matrix\n",
    "\n",
    "# construct 3x4 block matrix\n",
    "\n",
    "A = [B C ; \n",
    "    D E]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Column and row interpretation of a matrix.** An $m × n$ matrix $A$ can be interpreted as a collection of $n$ $m$-vectors (its columns) or a collection of $m$ row vectors (its rows). Julia distinguishes between a matrix (a two-dimensional array) and an array of vectors. An array (or a tuple) of column vectors can be converted into a matrix using the horizontal concatenation function `hcat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Array{Float64,1},1}:\n",
       " [1.0, 2.0]\n",
       " [4.0, 5.0]\n",
       " [7.0, 8.0]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [ [1., 2.], [4., 5.], [7., 8.] ] # array of 2-vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Float64,2}:\n",
       " 1.0  4.0  7.0\n",
       " 2.0  5.0  8.0"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = hcat(a...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `...` operator in `hcat(a...)` splits the array `a` into its elements, i.e., `hcat(a...)` is the same as `hcat(a[1], a[2], a[3])`, which concatenates `a[1], a[2], a[3]` horizontally. Similarly, `vcat` concatenates an array of arrays vertically. This is useful when constructing a matrix from its row vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Array{Float64,2},1}:\n",
       " [1.0 2.0]\n",
       " [4.0 5.0]\n",
       " [7.0 8.0]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [ [1. 2.], [4. 5.], [7. 8.] ] # array of 1x2 matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Float64,2}:\n",
       " 1.0  2.0\n",
       " 4.0  5.0\n",
       " 7.0  8.0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = vcat(a...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 Zero and identity matrices\n",
    "**Zero matrices.** A zero matrix of size m× n is created using zeros(m,n)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.0  0.0\n",
       " 0.0  0.0"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "zeros(2,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Identity matrices.** Identity matrices in Julia can be created many ways, for example by starting with a zero matrix and then setting the diagonal entries to one.\n",
    "\n",
    "The `LinearAlgebra` package also contains functions for creating a special identity matrix object `I`, which has some nice features. You can use `1.0*Matrix(I,n,n)` to create an `n× n` identity matrix. (Multiplying by $1.0$ converts the matrix into one with numerical entries; otherwise it has $Boolean$ entries.) This expression is pretty unwieldy, so we can define a function `eye(n)` to generate an identity matrix. This function is in the `VMLS` package, so you can use it once the package is installed. (The name `eye` to denote the identity matrix `I` traces back to the `MATLAB` language.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×4 Array{Float64,2}:\n",
       " 1.0  0.0  0.0  0.0\n",
       " 0.0  1.0  0.0  0.0\n",
       " 0.0  0.0  1.0  0.0\n",
       " 0.0  0.0  0.0  1.0"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eye(n) = 1.0*Matrix(I,n,n)\n",
    "eye(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Julia’s identity matrix `I` has some useful properties. For example, when it can deduce its dimensions, you don’t have to specify it. (This is the same as with common mathematical notation; see VMLS page [113](https://web.stanford.edu/~boyd/vmls/vmls.pdf#section*.135).)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Int64,2}:\n",
       " 1  -1   2\n",
       " 0   3  -1"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [ 1 -1 2; 0 3 -1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×5 Array{Int64,2}:\n",
       " 1  -1   2  1  0\n",
       " 0   3  -1  0  1"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[A I]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5×3 Array{Int64,2}:\n",
       " 1  -1   2\n",
       " 0   3  -1\n",
       " 1   0   0\n",
       " 0   1   0\n",
       " 0   0   1"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[A ; I]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Int64,2}:\n",
       " 1  2\n",
       " 3  4"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = [ 1 2 ; 3 4 ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Int64,2}:\n",
       " 2  2\n",
       " 3  5"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B + I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Ones matrix.** In VMLS we do not have notation for a matrix with all entries one. In Julia, such a matrix is given by `ones(m,n)`. \n",
    "\n",
    "**Diagonal matrices.** In standard mathematical notation, `diag(1, 2, 3)` is a diagonal $3 × 3$ matrix with diagonal entries $1, 2, 3$. In Julia such a matrix is created using the function `diagm`, provided in the `LinearAlgebra` package. To construct the diagonal matrix with diagonal entries in the vector `s`, you use `diagm(0 => s)`. This is fairly unwieldy, so the `VMLS` package defines a function `diagonal(s)`. (Note that you have to pass the diagonal entries as a vector.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "diagonal (generic function with 1 method)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    " diagonal(x) = diagm(0 => x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Int64,2}:\n",
       " 1  0  0\n",
       " 0  2  0\n",
       " 0  0  3"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diagonal([1,2,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A closely related Julia function `diag(X)` does the opposite: It takes the diagonal entries of the (possibly not square) matrix `X` and puts them into a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×4 Array{Int64,2}:\n",
       " 0   1  -2  1\n",
       " 2  -1   3  0"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = [0 1 -2 1; 2 -1 3 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       "  0\n",
       " -1"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diag(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Random matrices.** A random $m×n$ matrix with entries between $0$ and $1$ is created using `rand(m,n)`. For entries that have a normal distribution, `randn(m,n)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×3 Array{Float64,2}:\n",
       " 0.41485  0.482054  0.912613\n",
       " 0.49318  0.850034  0.910118"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rand(2,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Float64,2}:\n",
       "  0.311996  -0.342588\n",
       " -1.36241    0.453486\n",
       "  0.935138  -0.59365 "
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "randn(3,2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices. Functions for creating and manipulating sparse matrices are contained in the SparseArrays package, which must be installed; see page ix. Sparse matrices are stored in a special format that exploits the property that most of the elements are zero. The sparse function creates a sparse matrix from three arrays that specify the row indexes, column indexes, and values of the nonzero elements. The following code creates a sparse matrix\n",
    "\n",
    "$$\n",
    "  A = \n",
    "  \\begin{bmatrix} \n",
    "  −1.11 &  0.00 &  1.17 &  0.00 &  0.00 \\\\\n",
    "  0.15  & −0.10 &  0.00 &  0.00 &  0.00 \\\\\n",
    "  0.00  & 0.00  & −0.30 &  0.00 &  0.00 \\\\\n",
    "  0.00  & 0.00  & 0.00  & 0.13  & 0.00 \n",
    "  \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SparseArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 SparseMatrixCSC{Float64,Int64} with 6 stored entries:\n",
       "  [1, 1]  =  -1.11\n",
       "  [2, 1]  =  0.15\n",
       "  [2, 2]  =  -0.1\n",
       "  [1, 3]  =  1.17\n",
       "  [3, 3]  =  -0.3\n",
       "  [4, 4]  =  0.13"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = [ 1, 2, 2, 1, 3, 4 ] # row indexes of nonzeros\n",
    "N = [ 1, 1, 2, 3, 3, 4 ] # column indexes\n",
    "V = [ -1.11, 0.15, -0.10, 1.17, -0.30, 0.13 ] # values\n",
    "A = sparse(M, N, V, 4, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nnz(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices can be converted to regular non-sparse matrices using the `Array` function. Applying `sparse` to a full matrix gives the equivalent sparse matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 SparseMatrixCSC{Float64,Int64} with 4 stored entries:\n",
       "  [1, 1]  =  1.0\n",
       "  [3, 1]  =  2.0\n",
       "  [2, 2]  =  3.0\n",
       "  [1, 3]  =  4.0"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sparse([1, 3, 2, 1], [1, 1, 2, 3],\n",
    "    [1.0, 2.0, 3.0, 4.0], 3, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 1.0  0.0  4.0\n",
       " 0.0  3.0  0.0\n",
       " 2.0  0.0  0.0"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B = Array(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 SparseMatrixCSC{Float64,Int64} with 3 stored entries:\n",
       "  [1, 1]  =  1.0\n",
       "  [3, 1]  =  2.0\n",
       "  [2, 2]  =  3.0"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B[1,3] = 0.0;\n",
    "sparse(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sparse $m × n$ zero matrix is created with spzeros(m,n). To create a sparse $n × n$ identity matrix in Julia, use `sparse(1.0I,n,n)`. This is not a particularly natural syntax, so we define a function `speye(n)` in the `VMLS` package. The VMLS package also includes the function `speye(n)` which creates a sparse $n × n$ identity matrix, as well as `spdiagonal(a)`, which creates a sparse diagonal matrix with the entries of the vector a on its diagonal. \n",
    "\n",
    "A useful function for creating a random sparse matrix is sprand(m,n,d) (with entries between $0$ and $1$) and `sprandn(m,n,d)` (with entries that range over all numbers). The first two arguments give the dimensions of the matrix; the last one, `d`, gives the density of nonzero entries. The nonzero entries are chosen randomly, with about $mnd$ of them nonzero. The following code creates a random $10000 × 10000$ sparse matrix, with a density $10^{−7}$. This means that we’d expect there to be around $10$ nonzero entries. (So this is a *very* sparse matrix!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000×10000 SparseMatrixCSC{Float64,Int64} with 7 stored entries:\n",
       "  [2633 ,    90]  =  0.408501\n",
       "  [33   ,   723]  =  0.404502\n",
       "  [8559 ,   824]  =  0.496639\n",
       "  [5953 ,  3906]  =  0.982507\n",
       "  [7401 ,  3973]  =  0.631249\n",
       "  [3978 ,  4172]  =  0.753699\n",
       "  [3427 ,  8290]  =  0.248351"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sprand(10000,10000,10^-7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.3 Transpose, addition, and norm\n",
    "**Transpose.** In VMLS we denote the transpose of an $m × n$ matrix $A$ as $A^T$ . In Julia, the transpose of $A$ is given by `A'`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×4 Array{Int64,2}:\n",
       " 0   1  -2  1\n",
       " 2  -1   3  0"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = [0 1 -2 1; 2 -1 3 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×2 Adjoint{Int64,Array{Int64,2}}:\n",
       "  0   2\n",
       "  1  -1\n",
       " -2   3\n",
       "  1   0"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Addition, subtraction, and scalar multiplication.** In Julia, addition and subtraction of matrices, and scalar-matrix multiplication, both follow standard mathematical notation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Int64,2}:\n",
       " 0  4\n",
       " 7  0\n",
       " 3  1"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = [ 0 4; 7 0; 3 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Int64,2}:\n",
       " 1  2\n",
       " 2  3\n",
       " 0  4"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = [ 1 2; 2 3; 0 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Int64,2}:\n",
       " 1  6\n",
       " 9  3\n",
       " 3  5"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U+V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×2 Array{Float64,2}:\n",
       "  0.0  8.8\n",
       " 15.4  0.0\n",
       "  6.6  2.2"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.2*U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(We can also multiply a matrix on the right by a scalar.) \n",
    "\n",
    "Julia supports some operations that are not standard mathematical ones. For example, in Julia you can add or subtract a constant from a matrix, which carries out the operation on each entry. \n",
    "\n",
    "**Elementwise operations.** The syntax for elementwise vector operations described on page 10 carries over naturally to matrices. We add a period before a binary operator to change the interpretation to elementwise. For example, if $A$ and $B$ are matrices of the same size, then `C = A .* B` creates a matrix of the same size with elements $C_{ij} = A_{ij}B_{ij}$ . We can add a period after a function name to tell Julia that the function should be applied elementwise. Thus, if `X` is a matrix, then `Y = exp.(X)` creates a matrix of the same size, with elements $Y_{ij} = exp(X_{ij})$.\n",
    "\n",
    "**Matrix norm.** In $VMLS$ we use $‖A‖$ to denote the norm of an $m × n$ matrix, \n",
    "\n",
    "$$\n",
    "‖A‖ = \\left(\\sum_{i=1}^{m}\\sum_{j=1}^{n}A^{2}_{ij}\\right)^{\\frac{1}{2}}.\n",
    "$$\n",
    "\n",
    "In standard mathematical notation, this is more often written as $‖A‖F$ , where $F$ stands for the name Frobenius. In standard mathematical notation, $‖A‖$ usually refers to another norm of a matrix, that is beyond the scope of the topics in VMLS. In Julia, `norm(A)` gives the norm used in VMLS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5.5677643628300215, 5.5677643628300215)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [2 3 -1; 0 -1 4]\n",
    "norm(A),norm(A[:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Triangle inequality.** Let’s check that the triangle inequality $‖A+B‖ ≤ ‖A‖+‖B‖$ holds, for two specific matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.69041575982343, 7.795831523312719)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [-1 0; 2 2]; B= [3 1; -3 2];\n",
    "norm(A + B), norm(A) + norm(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.4 Matrix-vector multiplication\n",
    "In Julia, matrix-vector multiplication has the natural syntax `y=A*x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Int64,1}:\n",
       "  3\n",
       " -4"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [0 2 -1; -2 1 1]\n",
    "x = [2, 1, -1]\n",
    "A*x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Difference matrix.** An $(n − 1) × n$ difference matrix (equation ([6.5](https://web.stanford.edu/~boyd/vmls/vmls.pdf#equation.6.4.6)) of VMLS) can be constructed in several ways. A simple one is the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 Array{Float64,2}:\n",
       " -1.0   1.0   0.0  0.0\n",
       "  0.0  -1.0   1.0  0.0\n",
       "  0.0   0.0  -1.0  1.0"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "difference_matrix(n) = [-eye(n-1) zeros(n-1)] + [zeros(n-1) eye(n-1)];\n",
    "D = difference_matrix(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  1.0\n",
       "  2.0\n",
       " -1.0"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D*[-1,0,2,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since a difference matrix contains many zeros, this is a good opportunity to use\n",
    "sparse matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SparseArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×4 SparseMatrixCSC{Float64,Int64} with 6 stored entries:\n",
       "  [1, 1]  =  -1.0\n",
       "  [1, 2]  =  1.0\n",
       "  [2, 2]  =  -1.0\n",
       "  [2, 3]  =  1.0\n",
       "  [3, 3]  =  -1.0\n",
       "  [3, 4]  =  1.0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "difference_matrix(n) = [-speye(n-1) spzeros(n-1)] + [spzeros(n-1) speye(n-1)];\n",
    "D = difference_matrix(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "  1.0\n",
       "  2.0\n",
       " -1.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D*[-1,0,2,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Running sum matrix.** The running sum matrix (equation ([6.6](https://web.stanford.edu/~boyd/vmls/vmls.pdf#equation.6.4.6)) in VMLS) is a lower triangular matrix, with elements on and below the diagonal equal to one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "running_sum (generic function with 1 method)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function running_sum(n) # n x n running sum matrix\n",
    "    S = zeros(n,n)\n",
    "    for i=1:n\n",
    "        for j=1:i\n",
    "            S[i,j] = 1\n",
    "        end\n",
    "    end\n",
    "    return S\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " -1.0\n",
       "  0.0\n",
       "  2.0\n",
       "  2.0"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "running_sum(4)*[-1,1,2,0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative construction is `tril(ones(n,n))`. This uses the function `tril`, which sets the elements of a matrix above the diagonal to zero. Vandermonde matrix. An $m×n$ Vandermonde matrix (equation ([6.7](https://web.stanford.edu/~boyd/vmls/vmls.pdf#equation.6.4.7)) in VMLS) has entries $t^{j−1}_i$ for $i = 1, . . . ,m$ and $j = 1, . . . , n$. We define a function that takes an $m$-vector with elements $t_1, . . . , t_m$ and returns the corresponding $m × n$ Vandermonde matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vandermonde (generic function with 1 method)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function vandermonde(t,n)\n",
    "    m = length(t)\n",
    "    V = zeros(m,n)\n",
    "    for i=1:m\n",
    "        for j=1:n\n",
    "            V[i,j] = t[i]^(j-1)\n",
    "        end\n",
    "    end\n",
    "    return V\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 Array{Float64,2}:\n",
       " 1.0  -1.0  1.0   -1.0    1.0   \n",
       " 1.0   0.0  0.0    0.0    0.0   \n",
       " 1.0   0.5  0.25   0.125  0.0625\n",
       " 1.0   1.0  1.0    1.0    1.0   "
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vandermonde([-1,0,0.5,1],5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative shorter definition uses Julia’s `hcat` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "vandermonde (generic function with 1 method)"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vandermonde(t,n) = hcat( [t.^i for i = 0:n-1]... )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×5 Array{Float64,2}:\n",
       " 1.0  -1.0  1.0   -1.0    1.0   \n",
       " 1.0   0.0  0.0    0.0    0.0   \n",
       " 1.0   0.5  0.25   0.125  0.0625\n",
       " 1.0   1.0  1.0    1.0    1.0   "
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vandermonde([-1,0,0.5,1],5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.5 Complexity\n",
    "**Complexity of matrix-vector multiplication.** The complexity of multiplying an\n",
    "$m×n$ matrix by an $n$-vector is $2mn$ flops. This grows linearly with both $m$ and $n$.\n",
    "Let’s check this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.009904 seconds (7.85 k allocations: 450.466 KiB)\n",
      "  0.004259 seconds (5 allocations: 8.094 KiB)\n",
      "  0.043273 seconds (6 allocations: 39.297 KiB)\n",
      "  0.038998 seconds (6 allocations: 39.297 KiB)\n"
     ]
    }
   ],
   "source": [
    "A = rand(1000,10000); x = rand(10000);\n",
    "@time y = A*x;\n",
    "@time y = A*x;\n",
    "A = rand(5000,20000); x = rand(20000);\n",
    "@time y = A*x;\n",
    "@time y = A*x;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the second matrix-vector multiply, $m$ increases by a factor of $5$ and $n$ increases by a factor of $2$, so the complexity predicts that the computation time should be (approximately) increased by a factor of $10$. As we can see, it is increased by a factor around $7.4$. The increase in efficiency obtained by sparse matrix computations is seen from matrix-vector multiplications with the difference matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.049277 seconds (6 allocations: 78.359 KiB)\n",
      "  0.035362 seconds (51.89 k allocations: 2.720 MiB)\n"
     ]
    }
   ],
   "source": [
    "n = 10^4;\n",
    "D = [-eye(n-1) zeros(n-1)] + [zeros(n-1) eye(n-1)];\n",
    "x = randn(n);\n",
    "@time y=D*x;\n",
    "Ds = [-speye(n-1) spzeros(n-1)] + [spzeros(n-1) speye(n-1)];\n",
    "@time y=Ds*x;"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.1.1",
   "language": "julia",
   "name": "julia-1.1"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
